Εργαστηριακή Άσκηση 2

…

from scipy import signal
# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
# https://docs.scipy.org/doc/scipy/reference/signal.html
from scipy.fft import fft, fftfreq
from dash import Dash, dcc, html, Input, Output
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
import pandas
import plotly
import plotly.express as px
import plotly.graph_objs as go
# Δείτε την έκδοση της αριθμητικής βιβλιοθήκης numpy
import dash_bootstrap_components as dbc
import warnings
warnings.filterwarnings('ignore')
print("Libraries added successfully!")

Libraries added successfully!

import warnings
warnings.filterwarnings('ignore')

Μέρος2: Σχεδιασμός και υλοποίηση φίλτρων

Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων: α) τη μέθοδο των παραθύρων και β) τη μέθοδο των ισοϋψών κυματώσεων τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.

Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα).

Εδώ θα πειραματιστούμε με δύο σήματα: (i) το sonar του παραδείγματος, το οποίο εδώ διαβάζεται από ένα .txt αρχείο (έχει προέλθει με εξαγωγή του s από το MATLAB) και (ii) ένα σήμα μουσικής, το violin.wav (σήμα από μουσική βιολιού), το οποίο περιέχει υψηλότερες συχνότητες και έχει προέλθει με δειγματοληψία στα Fs_viol=44100 Hz.

Σήμα sonar

# Ανάγνωση δειγμάτων σήματος από txt file
with open('files/sima.txt') as f:
    s = [float(x) for x in f]
s=np.array(s)   
print('μέγεθος σήματος=', s.shape)
Fs=8192

μέγεθος σήματος= (6565,)

Στο πεδίο του χρόνου

…

t=np.arange(0,len(s))/Fs
# Time domain plot of x
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=s, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of x', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='t (sec)',
                  yaxis_title='Amplitude',
                  template='plotly_white',
                  xaxis=dict(range=[0, 0.8]),
                  yaxis=dict(range=[-0.05, 0.05]))
fig.show()

Ακούμε το σήμα

# Πρέπει να έχουμε εγκατεστημένη τη βιβλιοθήκη sounddevice
import sounddevice as sd
import ipywidgets as widgets
from IPython.display import display

def play_sound(b):
    sd.play(20*s,Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

Φάσμα (spectrum)

f, Pxx_den = signal.welch(s, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the PSD using a logarithmic scale on the y-axis
fig.add_trace(go.Scatter(x=f, y=Pxx_den, mode='lines',line=dict(color='#1F77B4')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα σήματος sonar',title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True),  # Customize grid color
)

# Show figure
fig.show()

Σήμα βιολιού

from scipy import signal
import scipy.io.wavfile
import numpy as np
import matplotlib.pyplot as plt

f=open("files/violin.wav", 'rb')
Fs_viol, s_viol = scipy.io.wavfile.read(f)
print('Fs_viol=',Fs_viol, ' number of samples=',len(s_viol))
f.close()

Fs_viol= 44100  number of samples= 220500

Στο πεδίο του χρόνου

…

tvl = np.arange(0, len(s_viol)) / Fs_viol

# Create Plotly figure
fig = go.Figure()

# Add trace for the signal
fig.add_trace(go.Scatter(x=tvl, y=s_viol, mode='lines',line=dict(color='#1F77B4')))

# Update layout
fig.update_layout(
    title='Signal over time',title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Time (s)',
    yaxis_title='Amplitude',
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show()

import ipywidgets as widgets
from IPython.display import display
import sounddevice as sd
import threading

# Define your sound and its sampling rate
# s_viol and Fs_viol should be defined before this point

def play_sound():
    sd.play(s_viol, Fs_viol)

def on_button_clicked(b):
    threading.Thread(target=play_sound).start()

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(on_button_clicked)

# Display the button
display(play_button)

Φάσμα (spectrum) και Φασματόγραμμα (spectorgram)

f, Pxx_den = signal.welch(s_viol, Fs_viol, nperseg=1024, noverlap=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the PSD using a logarithmic scale on the y-axis
fig.add_trace(go.Scatter(x=f, y=Pxx_den, mode='lines',line=dict(color='#1F77B4')))

# Update layout for a semilog plot
fig.update_layout(
    title='Power Spectral Density', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα [Hz]',
    yaxis_title='Πυκνότητα φάσματος ισχύος [V<sup>2</sup>/Hz]',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(),
    yaxis=dict(),  # Customize grid color
)

# Set y-axis range
fig.update_yaxes(range=[np.log10(0.5e-2), np.log10(1e5)])

# Show figure
fig.show()

f, tsp, Sxx = signal.spectrogram(s, Fs)

# Create Plotly figure
fig = go.Figure(data=go.Heatmap(
    z=10 * np.log10(Sxx),  # Convert power to dB
    x=tsp,  # Time for x-axis
    y=f,  # Frequency for y-axis
))

# Update layout
fig.update_layout(
    title='Spectrogram', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Χρόνος [sec]',
    yaxis_title='Συχνότητα [Hz]',
    template='plotly_white'
)

# Show figure
fig.show()

Βαθυπερατά φίλτρα

Η μέθοδος των παραθύρων

… …

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt


H = np.hstack((np.ones(int(Fs/8)), np.zeros(int(Fs-Fs/4)), np.ones(int(Fs/8))))

# Create the x values for the stem plot
x_values = np.arange(len(H))

# Create Plotly figure
fig = go.Figure()

# Add stems
for x, y in zip(x_values, H):
    fig.add_trace(go.Scatter(x=[x, x], y=[0, y], mode='lines', line=dict(color='#1F77B4')))

# Add markers
fig.add_trace(go.Scatter(x=x_values, y=H, mode='markers', marker=dict(color='#1F77B4')))

# Update layout
fig.update_layout(
    title='Stem plot', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Index',
    yaxis_title='Amplitude',
    template='plotly_white'
)

# Show figure
fig.show()


#its kinda slow

Ορθογωνικό παράθυρο (απλή περικοπή της h)

h=np.real(np.fft.ifft(H))
middle=int(len(h)/2)
h=np.hstack((h[middle:],h[:middle]))
h32=h[middle-16:middle+16]
h128=h[middle-64:middle+64]

x_values_h32 = np.arange(len(h32))

# Create Plotly figure for h32
fig_h32 = go.Figure()

# Add stems for h32
for x, y in zip(x_values_h32, h32):
    fig_h32.add_trace(go.Scatter(x=[x, x], y=[0, y], mode='lines', line=dict(color='#1F77B4')))

# Add markers for h32
fig_h32.add_trace(go.Scatter(x=x_values_h32, y=h32, mode='markers', marker=dict(color='#1F77B4')))

# Update layout for h32
fig_h32.update_layout(
    title='Stem plot of h32', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Index',
    yaxis_title='Amplitude',
    template='plotly_white'
)

# Show figure for h32
fig_h32.show()

freq,resp32 = signal.freqz(h32)
freq,resp128 = signal.freqz(h128)

# Create Plotly figure
fig = go.Figure()

# Add traces for h32 and h128 frequency responses
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp32), mode='lines', name='h32', line=dict(color='#1F77B4')))
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp128), mode='lines', name='h128', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας βαθυπερατού φίλτρου<br>(με ορθογωνικά παράθυρα)', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show()

Παράθυρα Hamming και Kaiser

h64=h[middle-32:middle+32]
freq,resp64 = signal.freqz(h64)
w_hamming=signal.hamming(len(h64))
h64_hamming = np.multiply(h64,w_hamming)
w_kaiser=signal.kaiser(len(h64),5)
h64_kaiser = np.multiply(h64,w_kaiser)

freq,resp64_hamming = signal.freqz(h64_hamming)
freq,resp64_kaiser = signal.freqz(h64_kaiser)

fig = go.Figure()

# Add trace for h64 (rectangular window)
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64), mode='lines', name='Rectangular', line=dict(color='#1F77B4')))

# Add trace for h64 with Hamming window
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64_hamming), mode='lines', name='Hamming', line=dict(color='green')))

# Add trace for h64 with Kaiser window
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64_kaiser), mode='lines', name='Kaiser', line=dict(color='red')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας βαθυπερατού φίλτρου<br>παράθυρα hamming (πράσινο) και kaiser (κόκκινο)<br>(το αντίστοιχο ορθογωνικό σε μπλε)', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show()

Φίλτρα ισοϋψών κυματώσεων

lpass = signal.remez(64, [0, 1000, 1300, Fs/2], [1, 0], fs=Fs)
freq,resp_pm = signal.freqz(lpass)

# Compute frequency response for the equiripple (Remez) filter
freq_pm, resp_pm = signal.freqz(lpass)

# Create Plotly figure
fig = go.Figure()

# Add trace for rectangular window response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp64), mode='lines', name='Rectangular Window', line=dict(color='#1F77B4')))

# Add trace for equiripple filter response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq_pm / np.pi, y=np.abs(resp_pm), mode='lines', name='Equiripple Filter', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας βαθυπερατού φίλτρου equirriple (πράσινο)<br>(το αντίστοιχο με ορθογωνικό παραθυρο σε μπλε)', title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show()

Εφαρμογή του φίλτρου

s_pm = signal.convolve(s,lpass,mode='same')/sum(lpass)

f, Pxx_den = signal.welch(s_pm, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the power spectral density of the filtered signal
fig.add_trace(go.Scatter(x=f, y=Pxx_den, mode='lines', line=dict(color='#1F77B4')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα φιλτραρισμένου σήματος sonar',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True, range=[np.log10(1e-15), np.log10(1e-7)])  # Adjust y-axis range
)

# Show figure
fig.show()

def play_sound(b):
    sd.play(s_pm, Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

Ζωνοπερατά φίλτρα

Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο

# Με αναλυτικό υπολογισμό της κρουστικής απόκρισης και παράθυρο kaiser
f1=800; f2=1600;  
Ts=1/Fs
f2m1=(f2-f1); f2p1=(f2+f1)/2; N=256
t=np.arange(-(N-1),N-1,2)*Ts/2
hbp=2/Fs*np.divide(np.multiply(np.cos(2*np.pi*f2p1*t),np.sin(np.pi*f2m1*t))/np.pi,t)
hbpw=np.multiply(hbp,signal.kaiser(len(hbp),5))

s_bp=signal.convolve(s,hbp,'same')

f, Pxx_den = signal.welch(s_bp, Fs, noverlap=128, nperseg=256)

# Create Plotly figure
fig = go.Figure()

# Add trace for the power spectral density
fig.add_trace(go.Scatter(x=f, y=Pxx_den, mode='lines', line=dict(color='blue')))

# Update layout for a semilog plot
fig.update_layout(
    title='Φάσμα ζωνοπερατού σήματος sonar',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Πυκνότητα φάσματος ισχύος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True, range=[np.log10(1e-15), np.log10(1e-7)])  # Adjust y-axis range
)

# Show figure
fig.show()

def play_sound(b):
    sd.play(20 * s_bp, Fs)

# Create a button widget
play_button = widgets.Button(description="Play Sound")

# Link the button to the function that plays the sound
play_button.on_click(play_sound)

# Display the button
display(play_button)

Ζωνοπερατό ισουψών κυματώσεων

bpass = signal.remez(128, [0, f1*0.9, f1*1.1, f2*0.95, f2*1.05, Fs/2], [0, 1, 0], fs=Fs)
freq,resp_pm = signal.freqz(bpass)

# Compute frequency response for the equiripple bandpass filter
freq, resp_pm = signal.freqz(bpass)

# Create Plotly figure
fig = go.Figure()

# Add trace for equiripple bandpass filter response
fig.add_trace(go.Scatter(x=0.5 * Fs * freq / np.pi, y=np.abs(resp_pm), mode='lines', name='Equiripple Bandpass', line=dict(color='green')))

# Update layout for a semilog plot
fig.update_layout(
    title='Απόκριση συχνότητας ζωνοπερατού φίλτρου equirriple',title_x=0.5,title_font=dict(size=20, color='black', family="Arial, sans-serif"),
    xaxis_title='Συχνότητα (Hz)',
    yaxis_title='Κέρδος',
    yaxis_type='log',  # Set y-axis to logarithmic scale
    template='plotly_white',
    xaxis=dict(showgrid=True),
    yaxis=dict(showgrid=True)
)

# Show figure
fig.show()